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DEVELOPMENT OF A BENCHMARK EXAMPLE FOR DELAMINATION FATIGUE 

GROWTH PREDICTION 


Ronald Krueger* 

ABSTRACT 

The development of a benchmark example for cyclic delamination growth prediction is 
presented and demonstrated for a commercial code. The example is based on a finite 
element model of a Double Cantilever Beam (DCB) specimen, which is independent of the 
analysis software used and allows the assessment of the delamination growth prediction 
capabilities in commercial finite element codes. First, the benchmark result was created for 
the specimen: The number of cycles to delamination onset was calculated from the material 
data for mode I fatigue delamination growth onset. Then, the number of cycles during stable 
delamination growth was obtained incrementally from the material data for mode I fatigue 
delamination propagation. For the benchmark case, where the results for delamination onset 
and growth were combined, the delamination length was calculated for an increasing total 
number of load cycles. Second, starting from an initially straight front, the delamination was 
allowed to grow under cyclic loading in a finite element model of a commercial code. The 
number of cycles to delamination onset and the number of cycles during stable delamination 
growth for each growth increment were obtained from the analysis. In general, good 
agreement between the results obtained from the growth analysis and the benchmark results 
could be achieved by selecting the appropriate input parameters. Overall, the results are 
encouraging but further assessment for mixed-mode delamination is required. 


1. INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common practice to 
characterize the onset and growth of delaminations. In order to predict delamination onset or 
growth, the calculated strain energy release rate components are compared to interlaminar fracture 
toughness properties measured over a range from pure mode I loading to pure mode II loading. 

The virtual crack closure technique (VCCT) is widely used for computing energy release rates 
based on results from continuum (2D) and solid (3D) finite element (FE) analyses and to supply the 
mode separation required when using the mixed-mode fracture criterion [1, 2]. The virtual crack 
closure technique was recently implemented into several commercial finite element codes. As new 
methods for analyzing composite delamination are incorporated into finite element codes, the need 
for comparison and benchmarking becomes important since each code requires specific input 
parameters unique to its implementation. 

An approach for assessing the delamination propagation capabilities in commercial finite 
element codes under static loading was recently presented and demonstrated for VCCT for 
ABAQUS® * 1 [3, 4], In these recent papers, benchmark results were created for full three-dimensional 
finite element models of the Double Cantilever Beam (DCB) and the Single Leg Bending (SLB) 
specimen. Then, starting from an initially straight front, the delamination was allowed to propagate 
in the finite element model. The load-displacement relationship and the total strain energy obtained 
from the propagation analysis results and the benchmark results were compared and good 


*R. Krueger, National Institute of Aerospace, 100 Exploration Way, Hampton, VA, 23666, resident at Durability, 
Damage Tolerance and Reliability Branch, MS 188E, NASA Langley Research Center, Hampton, VA, 23681, USA. 

1 ABAQUS 8 is a product of Dassault Systemes Simulia Corp. (DSS), Providence, RI, USA 


1 


agreements could be achieved by selecting the appropriate input parameters. Overall, the results 
were encouraging but it was determined that further assessment for mixed-mode delamination is 
required [3,4], 

The objective of the present study was to create a benchmark example, independent of the 
analysis software used, which allows the assessment of the delamination fatigue growth prediction 
capabilities in commercial finite element codes. At the beginning, a benchmark example is created 
based on incremental finite element models of a DCB specimen. To avoid unnecessary 
complications, experimental anomalies such as fiber bridging were not addressed. First, a sample 
material and cyclic loading was selected. Second, the number of cycles to delamination onset, No, 
was calculated from the mode I fatigue delamination growth onset data of the material. Third, the 
number of cycles during stable delamination growth, ANg, was obtained incrementally from the 
material data for mode I fatigue delamination propagation by using growth increments of Aa= 0.1 
mm. Fourth, the total number of growth cycles, Ng, was calculated by summing over the increments 
ANg • Fifth, the corresponding delamination length, a, was calculated by summing over the growth 
increments Aa. Finally, for the benchmark case where results for delamination onset and growth 
were combined, the delamination length, a, was calculated and plotted versus an increasing total 
number of load cycles Nt=Nd+Ng ■ After creating the benchmark, the approach was demonstrated 
for the commercial finite element code ABAQUS®. Starting from an initially straight front, the 
delamination was allowed to grow based on the algorithms implemented into the commercial finite 
element software. Input control parameters were varied to study the effect on the computed 
delamination increase during cyclic loading. It was assumed that delamination length increase 
during cyclic loading obtained from finite element analysis should closely match the growth shown 
in the benchmark example. The benchmark enabled the selection of the appropriate input 
parameters that yielded good agreement between the results obtained from the growth analysis and 
the benchmark results. Once the parameters have been identified, they may then be used with 
confidence to model delamination growth for more complex configurations. 


2. METHODOLOGY 

The methodology for delamination propagation, onset and growth was applied to the DCB 
specimen to create the benchmark example [5, 6], Since the required material property input data 
were not readily available in the open literature for a single material, a fictitious set of properties 
was constructed for this benchmarking exercise. Individual properties for commonly used 
graphite/epoxy tape materials were obtained from the open literature to create this set to represent a 
typical graphite epoxy composite. The material properties are given in Tables I and II. 

2.1 Static fracture toughness 

The mode I fracture toughness (mixed-mode ratio Gu/Gf=0 ) is generated experimentally using 
the Double Cantilever Beam (DCB) tests (as shown in Figure 1) [7], A fracture toughness G/ c = 
0.17 kJ/m 2 was used in this benchmarking exercise [8], 

2.2 Fatigue delamination growth onset 

The number of cycles to delamination onset, No, can be obtained from the delamination onset 
curve plotted in Figure 2 [9, 10]. The onset curve (solid green line) is a power law fit 

G = m 0 • N'p (1) 
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of the experimental data (open, green circles) obtained from a DCB test using the respective 
standard for delamination growth onset [9]. 


2.3 Fatigue delamination growth 

The number of cycles during stable delamination growth, Ng, can be obtained from the fatigue 
delamination propagation relationship (Paris Law) plotted in Figure 3 [8]. The delamination growth 
rate (solid purple line) can be expressed as a power law function 


da 

dN 


= c-G 1 


( 2 ) 


where da/dN is the increase in delamination length per cycle and G max is the maximum energy 
release rate at the front at peak loading. The factor c and exponent n are obtained by fitting the curve 
to the experimental data (open black circles) obtained from DCB tests [8]. The critical energy 
release rate or fracture toughness, G/ c , was included in the plot of Figure 3 (blue solid vertical line). 
Since composites do not exhibit the same threshold behavior commonly observed in metals, a cutoff 
value, G t h, was chosen below which delamination growth was assumed to stop (green solid vertical 
line) [8], It has to be noted that this benchmarking exercise ignores branching or fiber bridging and 
hence the Paris Law was not normalized with the static R-curve as recently suggested [11, 12], 


3. SPECIMEN AND FATIGUE TEST DESCRIPTION 

For the current numerical investigation, the Double Cantilever Beam (DCB) specimen, as 
shown in Figure 1, was chosen since it is simple, only exhibits the mode I opening fracture mode 
and had been used previously to develop an approach to assess the quasi-static delamination 
propagation simulation capabilities in commercial finite element codes [3, 4], To avoid unnecessary 
complications, experimental anomalies such as fiber bridging were not addressed. For the current 
study, a DCB specimen made of graphite/epoxy with an unidirectional layup, [0]24, was modeled. 
The material properties are given in Table I [8]. The material, layup, overall specimen dimensions 
including initial crack length, a, were identical to the specimen used earlier [3,4], 

For the cyclic loading of the specimen, guidance was taken from a draft standard designed to 
determine mode I fatigue delamination propagation [11], In the draft document, it is recommended 
to start the test at a maximum displacement, d max , which causes the energy release rate at the front, 
Gimax, to reach initially about 80% of G/ c 



= 0.8 


( 3 ) 


The maximum load, P ma x, and maximum displacement, cW/2, were calculated using the known 
quadratic relationship between energy release rate and applied load or displacement 

G P~ 

Imax 1 max 

c p 2 

' J '/c r crit 
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where P C ru and d cnt are the critical values. For the current study, a critical energy release rate 
G/ c =0. 17 kJ/m 2 was used and the critical values P cr it and 3 cr u (grey dashed lines) were obtained from 
the benchmark for static delamination propagation [3, 4] shown in the load-displacement plot in 
Figure 4. The calculated maximum load, P ma x, and calculated maximum displacement, 3 max /2, are 
shown in Figure 4 (dashed red line) in relationship to the static benchmark case (solid grey circles 
and dashed grey line) mentioned above. During constant amplitude cyclic loading of a DCB 
specimen under displacement control, the applied maximum displacement, 3 max /2=0.61 mm, is kept 
constant while the load drops as the delamination length increases (solid red circles and solid red 
line). The energy release rate corresponding to an applied maximum displacement 3 max l2=0.61 mm 
was calculated for different delamination lengths a using equation (5). The energy release rate 
decreases with increasing delamination length, a, as shown in Figure 5 (solid red circles and solid 
red line). Delamination growth was assumed to stop once the calculated energy release rate drops 
below the cutoff value, G t h, (green solid horizontal line). The static benchmark case (solid grey 
circles and dashed grey line in Figure 5), where the delamination propagates at constant Gi c (solid 
blue line in Figure 5) was included for comparison. 

In the ASTM draft document, it is suggested to use a load ratio 7?=0.1 for testing [11], The 
corresponding minimum load, P min , and minimum displacement, d min /2, were calculated 


P 

^ mm 


8 . 


= 0.1 


8 ,„,„ = 0 . 1 - 8 ,. 


( 6 ) 


Further, it was suggested in the ASTM draft document to use a frequency f= 10 Hz for testing [11]. 
A graphical representation of the cyclic fatigue loading is plotted in Figure 6. The applied 
displacement 8/2 is represented as a function of time, t 

d / 2 = [u 0 + b x • sin co(t - 1 0 )] ■ 3 max / 2 (7) 

where 3 max l2=0.61 mm is the maximum displacement. The constants ao= 0.55, bi= 0.45, the circular 
frequency co=20jt=62.832 and the starting time to= 0.025 are calculated from load ratio R=0. 1 and 
the frequency j= 10 Hz for testing. The resulting equation to calculate the applied displacement 3/2 
is shown in Figure 6. 


4. CREATING A BENCHMARK EXAMPLE FOR GROWTH PREDICTION 
4.1 Fatigue delamination growth onset 

The number of cycles to delamination onset, No, may be obtained by solving equation (1) for 

No. 


G = m 0 


■N 


in, 
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N 


D 



N D =c l -G c 
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where c, = — . Values for the constants c/ and C 2 are shown in Figure 7. 
m 1 

At the beginning of the test, the specimen is loaded initially so that the energy release rate at the 
front, Gimax, reaches about 80% of Gi c corresponding to Gi max = 0 . 1 362 kJ/rn 2 . The initial energy 
release rate is shown in the delamination onset plot of Figure 7 as a horizontal dashed red line. From 
the delamination onset curve, the number of cycles to delamination onset is determined as, /V/)=l 50 
shown as a vertical dashed red line. 

4.2 Fatigue delamination growth 

The number of cycles during stable delamination growth can be obtained by solving equation 
(2) for N (} 


N c = 



( 9 ) 


As mentioned above, the specimen is loaded initially so that the energy release rate at the front, 
Gimax, reaches about 80% of G/ c corresponding to Gy mav =0.1362 kJ/m 2 in the current study as shown 
in the Paris Law plot of Figure 8 (solid red square). 

For practical applications, equation (2) can be replaced by an incremental equivalent expression 


A a 
AN 


c-G'L 


( 10 ) 


where for the current study, increments of Aa= 0.1 mm were chosen. Starting at the initial 
delamination length a«= 30.5 mm, the energy release rates G l)tnca were obtained for each increment, i, 
from the curve fit (solid red circles and solid red line) plotted in Figure 5. These energy release rate 
values were then used to obtain the increase in delamination length per cycle or growth rate Act/ AN 
from the Paris Law in Figure 8. The growth rate rapidly decreases with increasing delamination 
length a as shown in Figure 9 (solid red circles and solid red line). The number of cycles during 
stable delamination growth, Ng, was calculated by summing the increments AN, 

No~t A^-i-GiL.-Aa (11) 

i-1 i=l C 


where k is the number of increments. The corresponding delamination length, a, was calculated by 
adding the incremental lengths Aa to the initial length ag. 

k 

a = a 0 + ^ Aa = a 0 + k ■ Act (12) 

1=1 


For stable delamination growth, the delamination length, a, 
increasing number of load cycles Ng (crosses and solid red line). 


is plotted in Figure 10 for an 
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4.3 Combined fatigue delamination onset and growth 

For the combined case of delamination onset and growth, the total life, Nt, may be expressed as 

N t =N d +N g (13) 

where, No, is the number of cycles to delamination onset and Ng, is the number of cycles during 
delamination growth [12], For this combined case, the delamination length, a, is plotted in Figure 1 1 
for an increasing number of load cycles Nt. For the first Nd cycles, the delamination length remains 
constant (horizontal red line), followed by a growth section where - over Ng cycles - the 
delamination length increases following the Paris Law (crosses and solid red line). Once a 
delamination length is reached where the energy release rate drops below the assumed cutoff value, 
G t h, (as shown in Figure 5) the delamination growth no longer follows the Paris Law (dashed grey 
line) and stops (horizontal solid red line). 

Applying the relationship 


t 



(14) 


where t is the time and / is the frequency, the development of the delamination length a can be 
plotted on a time scale assuming a frequency of 10 Hz as shown in Figure 12 (solid red line). 

4.4 Using the benchmark example to assess an automated analysis in a commercial FE code 

The load/displacement behavior of the DCB specimen - as shown in Figure 4 (solid red line) - 
can serve as an initial check for the finite element model. The correct input of model dimensions, 
material, layup and load application can thus be verified. During initial loading, the load and 
displacement should increase in a linear fashion and follow the initial slope until the maximum load, 
Pmax, and maximum displacement, d max , are reached. To minimize problems with numerical stability 
of the analysis, it is suggested that prescribed displacements, d max / 2 and d m J2 are applied in the 
analysis instead of nodal point loads, P max and P mm . The same approach was used to create the static 
benchmark case mentioned above [3,4]. Once the delamination starts to grow, the load is expected 
to drop while the applied maximum displacement cW/2=0.67 is expected to remain constant. 

The energy release rate may serve as an additional verification step before the delamination 
growth analysis is initiated. The computed values at the initial crack tip in a 2D model or along the 
delamination front in a 3D model should reach but not exceed the target value ofG/ maY = 0.136 kJ/m 2 . 
Once delamination growth starts in the model, the computed energy release rate should decrease 
with increasing delamination length, a, as shown in Figure 5. The curve fit (solid red line) can 
therefore be used to check the computed energy release rate during delamination growth. 

The growth rate Aa/AN, shown in Figure 9, decreases and the curve fit (solid red line) may be 
used as a check for the correct implementation of the Paris Law provided this output is available. 
For the delamination growth analysis, the delamination length, a, should increase with the number 
of cycles, N, as shown in Figure 10. The curve fit (solid red line) can therefore be used as a 
benchmark. 

A delamination length prediction analysis that accounts for delamination fatigue onset as well 
as stable growth should yield results that closely resemble the plot in Figure 11. The curve fit (solid 
red line) can therefore be used as a benchmark. 
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5. FINITE ELEMENT MODELING 

A typical two-dimensional finite element model of a Double Cantilever Beam (DCB) specimen 
is shown in Figure 13. The specimen was modeled with solid plane strain elements (CPE4) and 
solid plane stress elements (CPS4) in ABAQUS 8 Standard 6.8, 6.9 and 6.9EF. Along the length, all 
models were divided into different sections with different mesh refinement. The DCB specimen was 
modeled with six elements through the specimen thickness (2h) as shown in the detail of Figure 
13b. The resulting element length at the delamination tip was Aa=0.5 mm. A finer mesh, resulting in 
Aa=0.25 mm, was also generated as is shown in Figure 13c. Additionally, three coarser meshes with 
a reduced number of elements in the length direction were also generated as shown in Figures 14, 
resulting in Aa=1.0 mm (Figure 14b), Aa=1.25 mm (Figure 14c) and Aa=1.67 mm (Figure 14a). 

The plane of delamination was modeled as a discrete discontinuity in the center of the 
specimen. For the analysis with ABAQUS 8 6.8, 6.9 and 6.9EF, the models were created as separate 
meshes for the upper and lower part of the specimens with identical nodal point coordinates in the 
plane of delamination [13], Two surfaces (top and bottom surface) were defined to identify the 
contact area in the plane of delamination as shown in Figure 13b. Additionally, a node set was 
created to define the intact (bonded nodes) region. 

Typical three-dimensional finite element models of the DCB specimen are shown in Figures 15 
and 16. Along the length, all models were divided into different sections with different mesh 
refinement. A refined mesh was used in the center of the DCB specimen as shown in the detail of 
Figure 15b. Across the width, a uniform mesh was used to avoid potential problems at the transition 
between a coarse and finer mesh [3,4], Through the specimen thickness ( 2h ), six elements were 
used as shown in the detail of Figure 15b. The resulting element length at the delamination tip was 
Aa=0.5 mm. The specimen was modeled with solid brick elements (C3D8I) which had yielded 
excellent results in a previous studies [3,4]. Two coarser meshes with a reduced number of elements 
in the width and length directions, resulting in Aa=1.0 mm and Aa=2.0 mm, were also generated as 
shown in Figures 16a and b. 

Three models of the DCB specimen were generated with continuum shell elements (SC8R) as 
shown in Figures 17a to c. The continuum shell elements in ABAQUS 8 are used to model an entire 
three-dimensional body. Unlike conventional shells, which model a reference surface, the SC8R 
elements have displacement degrees of freedom only, use linear interpolation, and allow finite 
membrane defonnation and large rotations and, therefore, are suitable for geometric nonlinear 
analysis. The continuum shell elements are based on first-order layer-wise composite theory and 
include the effects of transverse shear deformation and thickness change [13]. In the x-y plane, the 
models had the same fidelity as the models made of solid brick elements C3D8I shown in Figures 
4b, 5a and 5b, resulting in an element length at the delamination tip Aa=0.5 mm, Aa=1.0 mm and 
Aa=2.0 mm, respectively. In the z-direction, only one element was used to model the thickness of 
the specimen. These less-refined models were used to study the effect on performance (CPU time), 
computed load/displacement behavior and growth prediction in comparison with the more refined 
models discussed above. 

For all the analyses perfonned, the low-cycle fatigue analysis in ABAQUS 8 Standard 6.8, 6.9 
and 6.9EF was used to model delamination growth at the interfaces in laminated composites [13, 
14]. A direct cyclic approach is part of the implementation and provides a computationally effective 
modeling technique to obtain the stabilized response of a structure subjected to constant amplitude 
cyclic loading. The theory and algorithm to obtain a stabilized response using the direct cyclic 
approach are described in detail in reference 14. Delamination onset and growth predictions are 
based on the calculation of the strain energy release rate at the delamination front using VCCT. To 
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determine propagation, computed energy release rates are compared to the input data for onset and 
growth from experiments as discussed in the methodology section. During the analysis, at least one 
element length at the crack tip is released along the interface after each stabilized cycle. 

For all analyses, the elastic constants, the input to define the fracture criterion, and the 
parameters for delamination onset and delamination growth (Paris Law) were kept constant. The 
elastic constants are given in Table I. The fracture toughness values and the parameters required for 
delamination onset and growth are given in Table II. The parameters to define the load frequency 
(/= 10 Hz), the load ratio (R=(). 1 ) as well as the minimum and maximum applied displacement 
(d m J 2=0.067 mm and <5 maI /2=0.67 mm) were also kept constant during all analyses. To study the 
effect on the computed onset and growth behavior during the analysis: 

• The number of terms used to define a Fourier series was varied. A Fourier series is 
used during the execution of the ABAQUS* Standard to approximate the periodic 
cyclic loading. 

• The size of the initial time increment used in the analysis was varied. 

• The input required to define the cyclic loading was altered. 

• The release tolerance was varied. Once a user specified release tolerance 
(( G-G c )/G c > release tolerance ) is exceeded during the analysis in ABAQUS® 
Standard, a cutback operation is performed which reduces the time increment. The 
cutback reduces the degree of overshoot and improves the accuracy of the local 
solution. 

• The solution controls were varied. 

It was assumed that the computed onset and growth behavior should closely match the 
benchmark results established earlier. Setting the value of the input parameters correctly is often 
an iterative procedure, which will be discussed later. Further details about the required input 
parameters are discussed in the appendix where sample input files are also provided. 


6. ANALYSIS 

6.1 Static analysis to verify correct input data and model response 

First, static analyses were preformed with propagation disabled to ensure that the input was 
correct and that the models responded as expected. For this static analysis, the applied displacement 
at the tip of the arms was ramped up to the maximum displacement d mar /2=0.67 for all models. To 
ensure that the model geometry and material input data for all models produced consistent results, 
the computed load-displacement behavior as shown in Figures 18 and 19 was evaluated. In 
comparison to the fatigue benchmark (solid grey circles and solid grey line) from Figure 4, which is 
based on analyses using solid 3D volume elements (C3D8I), the results are in good agreement. The 
model where plane strain elements (CPE4) were used exhibits a slightly more compliant behavior 
(dashed blue line) as shown in Figure 18. The model where plane stress elements (CPS4) were used 
exhibits a slightly stiffer behavior (dash-dot green line) whereas the results from the solid model 
(C3D8I) are as expected identical (solid black line) to the benchmark results. The results obtained 
from solid models with different mesh densities (see Figures 15 and 16) and from the continuum 
shell models (see Figure 17) are in good agreement with the benchmark as shown in Figure 19. 
Based on the results it was assumed that the geometry, maximum applied displacements and 
material input were defined correctly and all models adequately represented the benchmark case. 


In a second analysis step, only a single load cycle was analyzed, starting at the previously 
applied maximum displacement cW r /2=0.67. This step was performed to check that the amplitude 
input was defined correctly and resulted in the desired periodic cyclic loading during the analysis. 
Therefore, the displacement, d/2, - obtained as analysis output - was plotted versus the step time, t, 
as shown in Figure 20. The models where plane strain elements (CPE4) were used (blue crosses) 
yielded the same output as the models where plane stress elements (CPS4) were used (green x’s) 
and solid model (C3D8I) were used (open black circles). Based on a comparison with the desired 
fatigue loading (grey line) it is assumed that the amplitude input was defined correctly and the 
increments are small enough to adequately represent the desired periodic fatigue loading shown in 
Figure 6. 

Additionally, the computed mode I energy release rate - also obtained as analysis output - was 
plotted versus the step time, t, as shown in Figure 21. As desired, the energy release rate cycles 
between the expected maximum value Gi max = 0.136 kJ/m 2 and minimum value Gw=0.00136 kJ/m 2 . 
The model where plane strain elements (CPE4) were used yields lower energy release rates (blue 
crosses and solid blue line) while the model where plane stress elements (CPS4) were used yields 
higher energy release rates (green x’s and solid green line) which is consistent with the observations 
made above with respect to the slight variation in model stiffness. The results from the solid model 
(C3D8I), taken in the center of the specimen at y= 0.0, are as expected identical to the benchmark 
results (open black circles and solid black line). 

For the models made of 3D continuum elements (C3D8I and SC8R), the computed mode I 
strain energy release rate values were also computed for an applied maximum displacement 
6 max /2=0.67 and plotted versus the normalized width, y/B, of the specimen as shown in Figure 22. 
The results were obtained from models shown in Figures 15 through 17. Qualitatively, the mode I 
strain energy release rate is fairly constant in the center part of the specimen and drops progressively 
towards the edges. As desired, the energy release rate in the center of specimen reached the 
expected maximum value Gi max = 0.136 kJ/m 2 for all models used. Additionally, these results served 
as a verification that the VCCT procedure implemented in ABAQUS® 6.8 yielded the same results 
(open symbols) as an external post-processing routine (crosses) [4]. As expected, the mode II and 
mode III strain energy release rates were computed to be nearly zero and hence are not shown. 
Based on results, it was assumed that with respect to the computed energy release rates all models 
adequately represented the benchmark case. 

6.2 Results from fatigue onset and growth analysis 

In Figures 23 to 32, the delamination length, a, is plotted versus the number of cycles, N, for 
different input parameters and models. For all results shown, the analysis stopped when a 
10,000,000 cycles limit - used as input to terminate the analysis - was reached. 

6.2.1 Initial results 

Initial results, as plotted in Figures 23 and 24, were obtained using the specified default values 
as input parameters (see appendix). The results obtained from two-dimensional plane strain (open 
blue circles and solid blue line) and plane stress (open red squares and solid red line) models as well 
as full three-dimensional models (open green diamonds and solid green line) were within 4% of the 
benchmark results (grey crosses and solid grey line) as shown in Figure 23. For better visualization 
of the results and to be able to identify the differences in the results, the scale on the vertical axis 
was expanded as shown in Figure 24. For all results shown, the predicted onset occurs prior to the 
benchmark result onset, Ad= 150 cycles. The onset value is lowest for the plane stress results, 
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followed by the results obtained from 3D solid models and the plane strain results. This sequence 
can be explained by the computed energy release rate shown in Figure 21, where the values 
obtained from the plane stress model are slightly higher compared to the results obtained from 3D 
solid and plane strain models. During the growth phase of the analysis, the results lie on curves with 
nearly the same slope parallel to the benchmark which suggests that the Paris Law was implemented 
correctly and is - as expected - independent of the model. For all models, the threshold cutoff, where 
delamination growth is terminated and the delamination length remains constant, is predicted close 
to the number of cycles defined by the benchmark. 

6.2.2 Variation of input parameters 

Input parameters were varied to study the effect on the computed onset and growth behavior 
during the analysis and results are plotted in Figures 25 to 27. For this parametric study, only 
models with a refined mesh (Aa=0.25 mm, see Figure 13c) made of plane strain elements (CPE4) 
were used. Once a set of parameters was established that yielded good results, the effects of mesh 
size and element type on the results were studied. 

First, the effect of the initial time increment used in the analysis settings was studied as 
shown in Figure 25 (see appendix for details). The initial time increment was varied between 
z'o=0.01 (one tenth of a single loading cycle, 4=0. Is, open blue circles and solid blue line) and 
z'o=0.0001 (one thousandth of a single loading cycle, open green diamonds and solid green line). 
For larger initial time increments, the onset of delamination shifted towards a lower number of 
cycles. Reducing the initial time increments, however, significantly increased the computation 
time. Based on the results, it was therefore decided to use an initial time increment of z'o=0.001 
(open red squares and solid red line) for the remainder of the study to save computation time. This 
step is justified by the fact that the results obtained for z7f=0.001 were almost identical to the values 
obtained from the analysis where a smaller initial time increment was used (z'o=0.0001). 

Second, the input to define the cyclic loading was varied. In order to define the cyclic applied 
displacement, d/2, a set of parameters need to be determined as shown in the equation in Figures 26 
and 27 (see appendix for details). Selecting Ar=0 results in a sine curve representation of the cyclic 
load (analysis results are shown in Figure 26). Selecting Bj = 0 results in a cosine curve 
representation of the cyclic load (analysis results are shown in Figure 27). The selection of the 
starting time, to, causes a phase-shift. Further details about the input parameters are discussed in 
detail in the appendix where the corresponding cyclic applied displacements d/2 are plotted in 
Figure A8. As shown in Figures 26 and 27, the number of cycles to delamination onset is the most 
sensitive to input variations. The results obtained for the growth phase however, lie on curves with 
nearly the same slope, parallel to the benchmark. Also the threshold cutoff, where delamination 
growth is expected to stop and the delamination length remains constant, is predicted close to the 
number of cycles defined by the benchmark. A narrow band of results (green solid lines), was 
obtained when the number of terms used to define the Fourier series was increased to 50. A 
Fourier series is used in ABAQUS® Standard during the analysis to approximate the periodic 
cyclic loading (see appendix). The narrow band of results was in good agreement with the 
benchmark curve (grey solid line) compared to the results obtained for the default setting of 1 1 
Fourier terms which showed more scatter (black dashed lines). It should be noted that one term in 
a Fourier series is sufficient to exactly represent a simple periodic function such as the simple sine 
function used in the current example. At the same time, a large number of terms should improve the 
approximation of a more complicated cyclic load as demonstrated. It should also be noted that good 
results could only be obtained when the analyses were performed with ABAQUS 8 Standard 6.9EF. 
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Earlier versions showed unexplainable wide variations in results. Based on the current results, it was 
decided to use a sine curve representation (Aj= 0) in combination with the starting time, to= 0.0 for 
the remainder of the study. 

Third, the release tolerance was varied. Once a user specified release tolerance (( G-G c )/G c > 
release tolerance ) is exceeded during the analysis in ABAQUS® Standard, a cutback operation is 
performed which reduces the time increment The cutback reduces the degree of overshoot and 
improves the accuracy of the local solution [13]. In the current study, varying the input between 0.2 
(the default value) and 0.01 did not have any effect on the computed onset and growth behavior 
during the analysis. It is assumed that the release tolerance only affects static delamination 
propagation and not growth under cyclic loading studied here. 

Fourth, the solution controls were varied. It is generally not required nor recommended to 
modify the solution controls in ABAQUS 8 Standard. Since some of the solution controls, however, 
were modified in the DCB example problem provided by ABAQUS® [15], it was decided to briefly 
address the effect on the computed onset and growth behavior during the analysis. The first input 
parameter that defines the iteration number at which the periodicity condition is first imposed, 
was not modified and the default value was kept. A series of analyses were performed where the 
other four input parameters were varied between 100 and 10' 4 . Changing the input did not have any 
effect on the computed onset and growth behavior during the analysis, however it significantly 
influenced the computation time. The analysis required only about a tenth of the computation time 
when the second input parameter was modified from the default value (5T0' 3 ) to values larger than 
10' 1 . Changing the other remaining three parameters did not have any effect for this example. 
Further details about the input parameters are discussed in the appendix where also sample input 
files are provided. 

6.2.3 Variation of mesh size and element type 

The results obtained for models with different mesh sizes and different element types are shown 
in Figures 28 to 33. For models made of plane strain elements (CPE4), the element length, Aa, at the 
delamination tip was varied as shown in Figures 13 and 14. The results obtained from the respective 
models are plotted in Figures 28 and 29. Excellent agreement with the benchmark curve (grey 
crosses and grey solid line) could be achieved for element lengths up to Aa= 1.25 mm, as shown in 
Figure 28. For element length Aa=\.61 mm (purple triangles and solid purple line), the predicted 
onset occurs for a slightly higher number of cycles. The observed mismatch is largely due to the 
increased element length, which causes the first growth step to be larger. The following growth 
increments are relatively large due to the increased element length (Aa=\.61 mm). However, during 
the growth phase of the analysis, the results from all models lie on curves with the same slope 
parallel to the benchmark, which suggests that the Paris Law was implemented correctly and is, as 
expected independent, of the mesh. The analyses were repeated with the goal to eliminate the onset 
part of the analysis and focus on the stable growth section. Input parameters were chosen so that the 
onset part was extremely short and stable growth basically started immediately. Further details 
about the input parameters are discussed in detail in the appendix. The results are shown in 
Figure 29. As mentioned above, excellent agreement with the benchmark curve (grey crosses and 
grey solid line) could be achieved for all element lengths except for element length Aa=\.61 mm 
(purple triangles and solid purple line). For all models, the threshold cutoff, where delamination 
growth is expected to stop and the delamination length remains constant, is predicted close to the 
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number of cycles defined by the benchmark. The total computation time was between 70 s for the 
coarsest mesh and 1030 s for the finest mesh as shown in Figure 30 1 . 

The analyses were repeated for models made of plane stress elements (CPS4), where the 
element length, Aa, at the delamination tip was varied as before. The results obtained from the 
respective models are plotted in Figure 31. As before for the plane strain models, excellent 
agreement with the benchmark curve (grey crosses and grey solid line) could be achieved for all 
element lengths except for element length Aa=1.67 mm (purple triangles and solid purple line). The 
total computation time took between 80 s for the coarsest mesh and 1050 s for the finest mesh 1 . The 
computation times are included in Figure 30 for comparison. 

The results obtained for models made of 3D solid and continuum shell elements (shown in 
Figures 15 to 17) are shown in Figure 32. Good agreement with the benchmark curve (grey crosses 
and grey solid line) could be achieved for small element lengths (Aa= 0.5 mm). For larger element 
lengths, the results started to deviate from the benchmark. However, during the growth phase of the 
analysis, the results from all models lie on curves with the same slope parallel to the benchmark as 
observed before. Results obtained from the continuum shell element models (blue dashed lines) 
were almost identical when compared to results obtained from corresponding full 3D solid models 
(green solid line) with the same element length, Aa. The computational effort, however, was 
reduced by a factor of 2.5 to 2.9 for analyses performed using the continuum shell element models. 
For the models made of 3D solid elements, the total computation time was between 5800 s for the 
coarsest mesh and about 7 days for the finest mesh. For the corresponding models made of 
continuum shell elements, the total computation time was between 2000 s for the coarsest mesh and 
about 2.8 days for the finest mesh 1 . The computation times are also included in Figure 30 for 
comparison. 

For final comparison, finite element analyses were repeated with two-dimensional and three- 
dimensional models with the same element length, Aa= 0.5 mm, at the delamination tip. The results 
obtained from two-dimensional plane strain (open blue circles and solid blue line) and plane stress 
(open red squares and solid red line) models as well as full 3D solid models (open green diamonds 
and solid green line) and continuum shell elements (orange x’s and solid orange line) were within 
1% of the benchmark curve (grey crosses and solid grey line) as shown in Figure 33. The results 
obtained from the continuum shell element model (orange x’s and solid orange line) were almost 
identical compared to results obtained from the full 3D solid model (open green diamonds and solid 
green line) as mentioned above. The results obtained from plane stress models (open red squares 
and solid red line) are close to the results obtained from three-dimensional models. For all results 
shown, the predicted onset occurs first for the plane stress results and last for the plane strain results. 
Also, the delamination length obtained from plane strain models are slightly lower. This observation 
can be explained by looking at the computed energy release rate in Figure 21. The plane stress 
model yields a slightly higher energy release rate compared to the 3D solid and the plane strain 
model. Delamination onset would therefore occur first in plane stress models as indicated by the 
results plotted in Figure 33. During the stable growth phase, however, the results for all models lie 
on curves with nearly the same slope parallel to the benchmark as mentioned earlier. For all models, 
the threshold cutoff, where delamination growth is assumed to stop and the delamination length 
remains constant, is predicted close to the number of cycles defined by the benchmark. 


1 CPU time on Dual-Core AMD Opteron(tm) Processor 8220 SE 
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7. SUMMARY AND CONCLUSIONS 

The development of a benchmark example for cyclic delamination growth prediction is 
presented and demonstrated for the commercial finite element code ABAQUS® Standard. The 
example is based on a finite element model of a Double Cantilever Beam (DCB) specimen, which is 
independent of the analysis software used and allows the assessment of the delamination growth 
prediction capabilities in commercial finite element codes. First, the development of a benchmark 
example for delamination fatigue growth prediction was presented step by step. The number of 
cycles to delamination onset was calculated from the material data for mode I fatigue delamination 
growth onset. The number of cycles during stable delamination growth was obtained incrementally 
from the material data for mode I fatigue delamination propagation. For the combined benchmark 
case of delamination onset and growth, the delamination length was calculated for an increasing 
total number of load cycles. Second, starting from an initially straight front, the delamination was 
allowed to grow under cyclic loading. The number of cycles to delamination onset and the number 
of cycles during stable delamination growth for each growth increment were obtained from the 
analysis. 

The results showed the following: 

• In general, good agreement between the results obtained from the growth analysis and 
the benchmark results could be achieved by selecting the appropriate input parameters. 
However, selecting the appropriate input parameters was not straightforward and 
often required an iterative procedure. 

• The onset prediction appeared much more sensitive to the input parameters than the 
growth prediction. 

• Consistent results were obtained when input parameters were selected such that 50 
terms in the Fourier series were used during the execution of ABAQUS® Standard to 
approximate the periodic cyclic loading. 

• Good agreement between analysis results and the benchmark could be achieved when 
the initial time increment used in the analysis was about one tenth of a single loading 
cycle. 

• Best results were obtained when a sine curve representation of the cyclic applied 
displacement was selected in combination with the starting time, to=0.0. 

• The release tolerance did not have an effect on the analysis or the computed results. 

• Accurately computing the onset and growth required fine meshes with an element 
length at the tip Aa< 1.0 mm. 

• The solution controls in ABAQUS® Standard had to be modified in order to reduce 
computation time. Even with carefully selected input parameters, the analyses for three- 
dimensional models of a simple DCB specimen required days. 

• Although implemented in ABAQUS® Standard 6.8, version 6.9EF was required to 
obtain consistently good results. 

• Improvements are needed to make this analysis applicable to real case scenarios such as 
more complex specimens or structural components. 

Overall, the results are promising. In a real case scenario, however, where the results are 
unknown, obtaining the right solution will remain challenging. Further studies are required which 
should include the assessment of the propagation capabilities in more complex mixed-mode 
specimens and on a structural level. 
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Assessing the implementation in one particular finite element code illustrated the value of 
establishing benchmark solutions since each code requires specific input parameters unique to its 
implementation. Once the parameters have been identified, they may then be used with confidence 
to model delamination growth for more complex configurations. 
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TABLE I. MATERIAL PROPERTIES [8], 


Unidirectional Graphite/Epoxy Prepreg 

E\\ = 139.4 GPa 

£22= 10.16 GPa 

£33 = 10.16 GPa 

vi2 = 0.30 

vj3 = 0.30 

V23 = 0.436 

Gl2 = 4.6 GPa 

Gi3=4.6 GPa 

G23 = 3.54 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes the ply principal 
axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 denotes the 
direction transverse to the fiber in the plane of the lamina and index 33 the direction perpendicular to the plane of the 
lamina. 


TABLE II. FRACTURE PARAMETERS. 

Fracture Toughness Data [8] - Figure A1 
Gj c = 0.17 kJ/m 2 Gjj c = G/// c =0.49 kJ/m 2 rj= 1 .62 

Delamination Growth Onset Data [10] - Figures 2, A3, 

/w o=0.2023 »?j=-0. 078924 

Delamination Growth Rate Data (Paris Law) [8] - Figures 3, A5 
Gic = 0.17 kJ/m 2 G t h = 0.06 kJ/m 2 

c=2.44 1 0 6 «=10.61 
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APPENDIX 

Delamination fatigue growth analysis in ABAQUS® 

Delamination growth at the interfaces in laminated composites subjected to cyclic loadings 
can be simulated in ABAQUS 8 by specifying the low-cycle fatigue criterion propagation analysis 
using the direct cyclic approach [13]. The interface along which the delamination (or crack) 
propagates must be indicated in the model using a fracture criterion definition. The onset and 
growth of fatigue delamination at the interfaces are characterized by the relative fracture energy 
release rate. The fracture energy release rates at the crack tips in the interface elements are 
calculated based on the virtual crack closure technique (VCCT). The low-cycle fatigue analysis 
in ABAQUS 8 is a quasi-static analysis on a structure subjected to sub-critical cyclic loading. The 
low-cycle fatigue analysis in ABAQUS 8 uses the direct cyclic approach to obtain the stabilized 
cyclic response of the structure directly. The direct cyclic analysis uses a combination of Fourier 
series and time integration of the nonlinear material behavior to obtain the stabilized cyclic 
response of the structure iteratively and therefore avoids the numerical expense associated with a 
transient analysis. The direct cyclic analysis in ABAQUS 8 is therefore suited for very large 
problems in which many load cycles must be applied to obtain the stabilized response. The direct 
cyclic analysis in ABAQUS®, however, is limited to geometrically linear behavior and fixed 
contact conditions. The theory and algorithm to obtain a stabilized response using the direct 
cyclic approach are described in detail in the ABAQUS 8 Theory Manual [14], 

Required input for ABAQUS® 

The input required to perform a delamination onset and growth analysis in ABAQUS® 
Standard is discussed in the following paragraphs. It is assumed that the reader is familiar with 
ABAQUS 8 Standard and the syntax used in the input file (.inp). The focus is therefore on the 
specific input that relates to delamination propagation, low-cycle fatigue and the direct cyclic 
approach in ABAQUS 8 . Two example input files are given at the end of this appendix to provide 
an overview of an entire analysis and assist the readers in creating their own analyses. 

Input for delamination propagation 

The interface along which the delamination (or crack) propagates must be indicated in the 
model using a fracture criterion definition: 

*DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERION, TYPE=fatigue , MIXED MODE BEHAVIOR=BK , TOLERANCE=<tol> 

<cl>, <c2>, <c3>, <c4>, <rl>, <r2>,<GIc>,<GIIc>, 

<GIIIc>, <eta> 

where vcct_top and vcct_bot are interface surfaces as shown in Figure A1 and <tol> is the 
tolerance within which the crack propagation criterion must be satisfied. The input parameters 
for the fracture criterion are obtained from the static mixed-mode failure criterion, the 
delamination growth onset criterion and the growth rate shown in Figures A2 -A5. 

The critical energy release rates <gic>,<giic>,<giiic> and the curve fit parameter <eta> are 
obtained from the mixed-mode failure criterion as shown in Figure A2. A quasi static mixed-mode 
fracture criterion is determined for a material by plotting the interlaminar fracture toughness, G c , 
versus the mixed-mode ratio, Gu/Gt as shown in Figure A2. The fracture criterion is generated 
experimentally using pure Mode I (Gn/Gf=Q) Double Cantilever Beam (DCB) tests (as shown in 
Figure 1) [7], pure Mode II (Gi/Gi=l) End-Notched Flexure (ENF) test, and Mixed Mode Bending 
(MMB) tests of varying ratios of Gy and Gn • For the material used in this study, the experimental 
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data (open, blue circles) and mean values (fdled, blue circles) are shown in Figure A2 [8]. A 2D 
fracture criterion was suggested by Benzeggah and Kenane [16] using a simple mathematical 
relationship between G, and Gu/Gt 


. (Al) 

In this expression, G/ c and Gn c are the experimentally-determined fracture toughness data for mode I 
and II as shown in Figure A2. The factor r/ was determined by a curve fit using the Levenberg- 
Marquardt algorithm in the KaleidaGraph™ graphing and data analysis software [17]. 

The parameters <cl>,<c2> are obtained by solving the law for growth onset (shown in Figure 
A3) for the number of cycles N as shown in equation (8) and illustrated in Figure A4 (black curve). 
In the case where a law for growth onset is not available or immediate onset is desirable in the 
analysis, parameters such as suggested in the ABAQUS' example problem [15] may be chosen as 
shown in Figure A3 (red box and curve). The parameters <c3>,<c4> are obtained directly from the 
Paris Law as shown in Figure A5. The parameter <rl> is calculated from the energy release rate 
cutoff value, G t h, and the fracture toughness, G/ c , as shown in Figure A5. To calculate the 
parameter <r2> the user needs to define an energy release rate upper limit, G p i, above which the 
fatigue crack will grow at an accelerated rate as shown in Figure A5. For the current benchmark 
example, G,,i was chosen to be 90% of the fracture toughness. 

Input for cyclic loading 

Defining a low-cycle fatigue analysis using the direct cyclic approach in ABAQUS* 
Standard requires the definition of an amplitude curve which describes the relative load 
magnitude: 

* amplitude , name=test , DEFINITION=PERIODIC 
<n> , <omega> , <t 0 > , <A 0 > 

<A 1 >,<B 1 > 

where test is the label to be used to refer to the amplitude curve. The parameters defined in the 
first line are the number of terms in the Fourier series, <n>, the circular frequency, <omega>, in 
radians per time, the starting time, <t 0 >, and the constant term in the Fourier series, <a 0 >, as 
shown in Figure A6. The parameters defined in the second line are the first coefficient of the 
cosine terms, <a l >, and the first coefficient of the sine terms, <Bi>, also shown in Figure A6. 

The amplitude curve is then referenced in the definition of the cyclic loading. In the current 
example, prescribed displacements were used to simulate the cyclic opening of the arms of the DCB 
specimen: 

* BOUNDARY , AMPLITUDE=test 
LFRONTP, 1, 1, 0.67 

LFRONTM, 1, 1, -0.67 

where lfrontp and lfrontm are node sets located at the tip of each ann where the displacements 
are applied as shown in Figure Al. The factor 0.67 is used to multiply the relative magnitude 
defined by the amplitude curve (shown in Figure A6) and obtain the applied cyclic displacement, 
d/2, as shown in Figure A7. 

The direct cyclic approach in ABAQUS* Standard is used to obtain the stabilized cyclic 
response of a structure directly: 
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*direct cyclic , fatigue 

< io > / <t S > / I / ^ n i' > / <n max > / ^An> , <i max >, 

/ / <N t > / / 

where the parameter fatigue is used to perform a low-cycle fatigue analysis. The parameters 
defined in the first line are the initial time increment, <i 0 >, the time of a single loading cycle, 
<t s > (as shown in Figure A7), the minimum time increment allowed (not used), the maximum 
time increment allowed (not used), the initial number of terms in the Fourier series, cn^, the 
maximum number of terms in the Fourier series, <n max >, the increment in number of terms in the 
Fourier series, <An> , and the maximum number of iterations allowed in a step, <i max >. The 
parameters defined in the second line are the minimum increment in number of cycles over 
which the damage is extrapolated forward (default used), the maximum increment in number of 
cycles over which the damage is extrapolated forward (default used), the total number of cycles 
allowed in a step, <n t >, and the damage extrapolation tolerance (default used). The time of a 
single loading cycle was kept constant at 4=0.1 s for all analyses. Most analyses were run up to 
Nf= 1 0 7 cycles in order to reach the threshold after which delamination growth stops as shown in 
Figures 1 1 and 12. All other input parameters were varied and the effect on the results studied as 
discussed in the main part of this report. 

Control parameters direct cyclic analysis 

Solution controls in ABAQUS® Standard can be reset and modified by using keyword 

*controls , type=direct cyclic 

<I PI >, <CR n > , <CU n > , <CR 0 > , <CU 0 > , 

where the parameter direct cyclic is used to set parameters that will be used to control the 
stabilized state and plastic ratcheting detections and to specify when to impose the periodicity 
condition for direct cyclic analysis. If this keyword is omitted, default parameters are used. The 
parameters defined in the first line are the iteration number at which the periodicity condition is 
first imposed, <i PI >, (default used), the stabilized state detection criterion for the ratio of the 
largest residual coefficient on any terms in the Fourier series to the corresponding average flux 
norm, <cr„>, the stabilized state detection criterion for the ratio of the largest correction to the 
displacement coefficient on any terms in the Fourier series to the largest displacement 
coefficient, <cu n >, plastic ratchetting detection criterion for the ratio of the largest residual 
coefficient on the constant term in the Fourier series to the corresponding average flux norm, 
<cr 0 >, and the plastic ratchetting detection criterion for the ratio of the largest correction to the 
displacement coefficient on the constant term in the Fourier series to the largest displacement 
coefficient, <cu 0 > [13]. These control parameters were varied and the effect on the results studied 
as discussed in the main part of this report. 


Example input files 

Two example input files are given to provide an overview of an entire analysis and assist the 
readers in creating their own analyses. The analysis was divided into two steps. In the first step, a 
small static preload step was introduced as a work around to avoid problems discovered with the 
initial contact conditions (ABAQUS R bug v68_1987). The second step was set-up to perform the 
desired cyclic analysis. It was found that the prescribed displacements in the static preload step 
had to be small (0.0067 mm) compared to the prescribed displacements (d„„„/2=0.067, 
(5 max /2=0.67) which were used to simulate the cyclic opening of the arms of the DCB specimen. 
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Analyses where larger preload steps were chosen did not converge in the second step and 
terminated prematurely. 

For all analyses, the input to define the fracture criterion (<Gic>,<Giic>,<Giiic>,<eta>), the 
parameters for de lamination onset (<cl>,<c2>), and delamination growth (Paris Law) 
(<c3>,<c4>, <rl>, <r2>) were kept constant. The parameters to define the load frequency 
(<omega> , <a 0 >) as well as the minimum and maximum applied displacement (<t s >, <5 max /2=0.67) 
were also kept constant during all analyses. Other parameters required to define the Fourier 
series which is used to define the cyclic load, the initial time increment and the optional input 
parameters to control the solution in ABAQUS" Standard were varied. The ABAQUS" keywords 
shown in bold type were discussed in detail in the previous paragraphs. 


Input file for fatigue onset and growth analysis 

*HEADING 

DCB-UD-T300/ 1076 , a=30.5 nun 
units: nun, N, MPa 

*** elements , nodes, material , etc 


*NSET , NSET=BONDED , GENERATE 
253, 1693, 8 

**** surface and contact definition for VCCT **** 

* SURFACE, TYPE=ELEMENT , NAME=VCCT BOT 
EL_BOT , S3 

* SURFACE, TYPE=ELEMENT , NAME=VCCTJTOP 
ELTOP, SI 

* CONTACT PAIR, INTERACTION=VCCT , AD JUST=BONDED , small sliding 
VCCTTOP , VCCTBOT 
* SURFACE INTERACTION, NAME=VCCT 
<width> 

* INITIAL CONDITION, TYPE=CONTACT 
VCCTTOP, VCCTBOT, BONDED 
* * 

*NSET , NSET=LFRONTP 
8 

*NSET , NSET=LFRONTM 
1 

**** VCCT fatigue input 
‘parameter 

** Damage and tolerance parameters 
tol=0 . 001 

** Fracture toughness : 

GIc = 0.17030 
GIIc = 0.49360 
GIIIc = 0.5 
** B-K parameter : 
eta=l . 62 

** width in the plane stress /strain direction 
width =25.0 

** fatigue crack growth data ** 
cl=2 . 8E-09 
c2=-12 . 415 
c3=2 . 44E+06 
c4=10 . 61 
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* * * * 


rl=0 . 353 
r2=0 . 9 

*** amplitude 
* amplitude , name=test , DEFINITION=PERIODIC 
1,62.832,0. ,0.55 
0,0.45 

*** history data *** 

*** static ramp up step **** 

*STEP, NLGEOM, INC= 10000 

*STATIC 
0.001, 0.001 
* * 

*DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 
* FRACTURE CRITERION, TYPE=VCCT , MIXED MODE BEHAVIOR=BK 
1 . 0e6 , 1 . 0e6 , 1 . 0e6 , <eta> 

* * 

* BOUNDARY, TYPE=DISPLACEMENT 
LFRONTP, 1, 1, 0.00067 
LFRONTM, 1, 1, -0.00067 

** field and history output ** 

*OUTPUT , FIELD, VARIABLE=PRESELECT , FREQ=1 
*Output , history , VARIABLE=PRESELECT , f req=l 
*NODE output, NSET=LFRONTP 
RF1 

*NODE output, NSET=LFRONTP 
U1 

*END STEP 
*STEP , INC= 10000 
*direct cyclic , fatigue 
0.001,0.1, , ,11,11,5,10, 

, ,10000000, , 

*DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERION, TYPE=fatigue , MIXED MODE BEHAVIOR=BK , TOLERANCE=<tol> 
<cl>, <c2>, <c3>, <c4>, <rl>, <r2>, <GIc>, <GIIc>, 

<GIIIc>, <eta> 

*** run analysis first with default values 
*controls , type=direct cyclic 
,100,5.E-3,5.E-3,5.E-3 
* BOUNDARY , AMPLITUDE=test 
LFRONTP, 1, 1, 0.67 

LFRONTM, 1, 1, -0.67 
* * 

** field and history output ** 

*OUTPUT , FIELD, VARIABLE= 

*ELEMENT OUTPUT 
cycleini , status , sdeg 

*CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
dbt , dbsf , dbs , openbc , crsts , enrrt , ef enrrtr , bdstat 
*0utput , history , VARIABLE=PRESELECT , f req=25 
*N0DE output, NSET=LFRONTP 
RF1 

*N0DE output, NSET=LFRONTP 
U1 

*END STEP 
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Input file to verify correct input data and model response 

It is recommended to perform a static analysis first, to verify the correct input parameters. A 
single cycle is performed by replacing the input for direct cyclic analysis 

*direct cyclic , fatigue 
0.001,0.1, , ,11,11,5,10, 

, , 10000000 , , 

in the second step with 

*STATIC 

0.001, 0.1, l.E-10, 0.001 

By monitoring the applied displacements, d/2, during the simulation, the input data for the cyclic 
analysis can be checked as shown in the example in Figure A8. For all analysis shown, the first step 
consisted of a static preload step up to 5/2=0.00067. In the second step, the number of terms in the 
Fourier series, <n>, the circular frequency, <omega>, the constant term in the Fourier series, <a 0 >, 
and the time of a single loading cycle, <t s > were kept the same for all analyses. Further, the factor 
0.67 used to multiply the relative magnitude defined by the amplitude curve to obtain the applied 
displacement, 5/2, was kept constant. The applied displacement, 5/2, therefore varied between 
dmaxl 2=0.67 mm and 8 m J 2=0.067 mm at a frequency of 10 Hz as desired. As an example, the 
starting time, <t 0 >, the first coefficient of the cosine terms, <a 2 >, and the first coefficient of the 
sine terms, <b 2 >, were varied to create different sine waves (in red) and cosine waves (in blue) 
with the same displacement maximum, minimum and frequency as shown in Figure A8. 


*HEADING 

DCB-UD-T300/ 1076 , a=30.5 mm 
units: mm, N, MPa 

*** elements , nodes, material , etc 


*NSET , NSET=BONDED , GENERATE 
253, 1693, 8 

**** surface and contact definition for VCCT **** 

* SURFACE, TYPE=ELEMENT , NAME=VCCT_BOT 
EL_BOT , S3 

* SURFACE, TYPE=ELEMENT , NAME=VCCTJTOP 
EL TOP , SI 

* CONTACT PAIR, INTERACTION=VCCT , AD JUST=BONDED , small sliding 
VCCTTOP , VCCTBOT 
* SURFACE INTERACTION, NAME = VCCT 
<width> 

* INITIAL CONDITION, TYPE=CONTACT 
VCCTTOP, VCCTBOT, BONDED 
* * ** 

*NSET , NSET=LFRONTP 
8 

*NSET , NSET=LFRONTM 
1 

**** VCCT fatigue input 
*parameter 

** Damage and tolerance parameters 
tol=0 . 001 

** Fracture toughness : 
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GIc = 0.17030 
GIIc = 0.49360 
GIIIc = 0.5 
** B-K parameter : 
eta=l . 62 

** width in the plane stress/strain direction 
width =25.0 

** fatigue crack growth data ** 
cl=2 . 8E-09 
c2=-12 . 415 
c3=2 . 44E+06 
c4=10 . 61 
rl=0 . 353 
r2=0 . 9 
* * * 

*** amplitude **** 

* amplitude , name=test , DEFINITION=PERIODIC 

1,62.832,0. ,0.55 

0,0.45 

*** history data *** 

*** static ramp up step **** 

*STEP , NLGEOM, INC= 10000 
* STATIC 
0.001, 0.1 
* * 

* DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 
* FRACTURE CRITERION, TYPE=VCCT , MIXED MODE BEHAVIOR=BK 
1 . 0e6 , 1 . 0e6 , 1 . 0e6 , <eta> 

* * 

* BOUNDARY, TYPE=DISPLACEMENT 
LFRONTP, 1, 1, 0.067 
LFRONTM, 1, 1, -0.067 

** field and history output ** 

*NODE PRINT, NSET=BONDED, GLOBAL=YES , FREQ=1 
COORD 

*OUTPUT , FIELD, VARIABLE=PRESELECT , FREQ=1 
*Output , history , VARIABLE=PRESELECT , f req=l 
*NODE output, NSET=LFRONTP 
RF1 

*NODE output, NSET=LFRONTP 
U1 

*CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP , NSET=BONDED 
dbt , enrrt , bdstat 

*END STEP 

*STEP, NLGEOM, INC= 10000 

*** run static analysis with one cycle first to verify input 
*STATIC 

0.001, 0.1, l.E-10, 0.001 

*DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 
* FRACTURE CRITERION, TYPE=VCCT , MIXED MODE BEHAVIOR=BK 
1 . 0e6 , 1 . 0e6 , 1 . 0e6 ,<eta> 

* BOUNDARY , AMPLITUDE=test 
LFRONTP, 1, 1, 0.67 
LFRONTM, 1, 1, -0.67 
* * 

** field and history output ** 

*OUTPUT , FIELD, VARIABLE=PRESELECT , FREQ=1 
*ELEMENT OUTPUT 
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status 

* CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
dbt , dbsf , dbs , openbc , crsts , enrrt , ef enrrtr , bdstat 

*Output , history , VARIABLE=PRESELECT , f req=l 

*NODE output, NSET=LFRONTP 
RF1 

*NODE output, NSET=LFRONTP 
U1 

* CONTACT OUTPUT, MASTER=VCCT BOT , SLAVE=VCCT_TOP , NSET=BONDED 
dbt , enrrt , bdstat 

*END STEP 
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dimensions 



B 25.0 mm 

2h 3.0 mm 

2L 150.0 mm 

a 0 30.5 mm 


fatigue loading 

W 2 067 mm 

6 min /2 0.067 mm 
R 0.1 
f 10.0 Hz 

layup: [0] 24 


Figure 1. Double Cantilever Beam Specimen (DCB) 
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G, kJ/m 2 


Figure 3. Delamination growth rate ( Paris Law) [8]. 
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Figure 4. Critical load-displacement behavior for a DCB specimen. 
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Figure 6. Cyclic fatigue loading for DCB specimen. 
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Figure 7. Delamination growth onset for DCB specimen. 
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Figure 8. Delamination growth rate ( Paris Law). 
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delamination length a, mm 

Figure 9. Delamination growth rate behavior for DCB specimen. 
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Figure 10. Stable delamination growth behavior for DCB specimen. 
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Figure 11. Delamination onset and growth behavior for DCB specimen. 
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a. Deformed FE-model ofDCB specimen with initial delamination before growth 



b. Detail of a FE-model of a DCB specimen 



c. Detail of a refined FE-model of a DCB specimen 
Figure 13. Two-dimensional finite element model of a DCB specimen. 
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a. Detail of a FE-model of a DCB specimen 


t 



b. Detail of a FE-model of a DCB specimen 



x 


c. Detail of a refined FE-model of a DCB specimen 
Figure 14. Details of different two-dimensional finite element models of a DCB specimen. 
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a. Deformed three-dimensional FE-model with initial delamination front before growth 



b. Detail of three-dimensional FE-model around delamination front 
Figure 15. Full three-dimensional finite element model of a DCB specimen. 
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a. Deformed three-dimensional model ofDCB specimen with a coarse mesh 



b. Deformed three-dimensional model of a DCB specimen with a coarse mesh 
Figure 16. Coarse full three-dimensional finite element models of a DCB specimen. 
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bonded 

nodes 


Aa=0.5 mm 

initial straight delamination front 


a. Deformed continuum-shell model ofDCB specimen with a fine mesh 



b. Deformed continuum- shell model of a DCB specimen with a coarse mesh 



c. Deformed continuum- shell model of a DCB specimen with a coarse mesh 
Figure 17. Continuum-shell finite element models of a DCB specimen. 
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load P, N 



applied opening displacement 5/2, mm 


Figure 18. Critical load-displacement behavior for a DCB specimen. 
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Figure 19. Critical load-displacement behavior for a DCB specimen. 
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Figure 21. Computed energy release rate. 
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Figure 22. Computed strain energy release rate distribution across the width of a DCB specimen. 
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Figure 23. Computed delamination onset and growth: Initial results. 
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Figure 24. Computed delamination onset and growth: Detail of initial results. 
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Figure 25. Computed delamination onset and growth obtained for different initial time increments. 
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delamination 
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Figure 26. Computed delamination onset and growth obtained for 
sine representation of the cyclic loading. 
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Figure 27. Computed delamination onset and growth obtained for 
cosine representation of the cyclic loading. 
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Figure 28. Computed delamination onset and growth behavior for different plane strain models. 
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Figure 29. Computed stable delamination growth behavior for different plane strain models. 
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crack tip element length Aa, mm 

Figure 30. Required analysis time for models with different crack tip element lengths. 
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Figure 31. Computed delamination onset and growth behavior for different plane stress models. 
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Figure 32. Computed delamination onset and growth behavior for different continuum models. 
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Figure 33. Computed delamination onset and growth behavior for different 2D and 3D models. 
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Figure A2. Mixed mode fracture criterion . 
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N d , cycle 

Figure A3. Delamination growth onset for DCB specimen. 


N, cycles 


10 7 
10 5 
10 3 
10 1 
10' 1 
10' 3 

0.00 0.05 0.10 0.15 0.20 0.25 



! *FRACTURE CRITERION, TYPE=f atigue, MIXED MODE BEHAVIOR=BK , 

: TOLERANCE=<tol> 

r <cl>,<c2>, <c3> , <c4> , <r 1> , <r2> , <GIc> , <GIIc> , 

: <GIIIc> , <eta> 

I I I I I I 1 1 I I I I I I I I I I I I I I I I 


G , kJ/m 2 (N/mm) 

max x ' 

Figure A4. Delamination growth onset input data. 
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Figure A5. Delamination growth rate (Paris Law) for T300/914C. 



Figure A6. Amplitude curve (relative load magnitude) . 
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Figure A7. Cyclic fatigue loading input. 
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Figure A8. Cyclic fatigue loading for DCB specimen. 
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